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Abstract 

Plasma instabilities can play a fundamental role in plasma equilibration. There are similarities 
and differences between plasma instabilities in abelian and non-abelian gauge theories. In par- 
ticular, it has been an open question whether non-abelian self-interactions are the limiting factor 
in the growth of non-abelian plasma instabilities. We study this problem with 3+1 dimensional 
numerical simulations. We find that non-abelian plasma instabilities behave very differently from 
abelian ones once they grow to be non-perturbatively large, in contrast with earlier results of 1+1 
dimensional simulations. In particular, they grow more slowly at late times, with linear rather 
than exponential dependence of magnetic energy on time. 



I. INTRODUCTION 



A fundamental problem in the theoretical study of quark-gluon plasmas (QGPs) is to un- 
derstand how such plasmas equilibrate. By what processes would a heavy ion collision first 
produce a quark-gluon plasma that is in approximate local equilibrium and/or expanding 
hydrodynamically? 1 At the very least, it would be useful and interesting to answer this ques- 
tion in a simplifying theoretical limit: that of arbitrarily high-energy collisions, where the 
running strong coupling a s at relevant scales is arbitrarily small due to asymptotic freedom. 
Even this weak-coupling limit is rich and complicated. Various authors have used weak- 
coupling techniques to study the initial creation and interaction of the non-perturbatively 
dense small- a; glue which will eventually develop into the quark-gluon plasma [2-8]. Baier 
et al. [9] have investigated the effects of two-particle scattering and Bremsstrahlung on sub- 
sequently equilibrating the plasma, once the initial small- x glue has expanded and diluted 
to perturbative densities, where it can be treated as a collection of individual gluons. How- 
ever, it now seems clear that the process of thermalization is not controlled solely by such 
individual particle collisions in the weak-coupling limit [10]; one must account for collective 
processes in the plasma in the form of plasma instabilities. The instabilities are known 
as Weibel or filamentary instabilities [11], and they have a very long history in traditional 
plasma physics. Their possible relevance to QGP equilibration has been proposed for roughly 
twenty years by Mrowczyriski and others [12-23]. 

In the case of electromagnetism, a simple example of a Weibel instability appears for 
two uniform, interpenetrating beams of charged particles traveling in opposite directions, as 
depicted in Fig. la. Crudely think of each beam as a superposition of many wires carrying 
current, and recall that parallel wires magnetically attract if their currents are aligned, or 
repel if opposite. The "wires" are therefore unstable to clumping, as depicted in Fig. lb. 
More thorough discussions of the qualitative origin of Weibel instabilities may be found in 
Refs. [20, 24] and [10]. Very roughly speaking, Weibel instabilities occur in collisionless 
plasmas whenever the velocity distribution of the plasma is anisotropic in its rest frame. 
(See Ref. [10] for a more precise statement.) For the purpose of analyzing the instability, 
the plasma may be regarded as collisionless whenever the distance and time scales associated 
with the instability are found to be small compared to those for individual, random collisions 
of particles in the plasma. 2 



(a) (b) 

FIG. 1: (a) Two uniform, inter-penetrating beams of charged particles, moving left and right; (b) 
the filamentation of those beams by the Weibel instability. 

There are magnetic fields associated with the filaments in Fig. lb. As the instability 



1 Hydrodynamic behavior does not necessarily require local thermal equilibrium, as discussed in Ref. [1]. 

2 See Ref. [10] for an analysis of this point in the context of the bottom-up scenario [9] for thermalization. 
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progresses, these fields become stronger and stronger. Growth of the instabilities is initially 
exponential in time. In electromagnetic plasmas, growth stops only when the magnetic fields 
become so large that they have a non-perturbatively large effect on the particle trajecto- 
ries. These effects of these large fields can drive isotropization of plasmas that are initially 
anisotropic [25]. This isotropization due to collective plasma phenomena is the basis of our 
own recent scenario (with Lenaghan) for the early onset of hydrodynamic behavior in the 
weak coupling limit [1], as compared to estimates based on individual 2-particle collisions. 
To analyze correctly the effects of plasma instabilities, however, it is crucial to understand 
whether the (chromo-) magnetic fields grow as large in non-abelian gauge theories as in 
abelian ones. Because non-abelian fields interact with each other, there are two possibilities 
for what might limit their growth. They might stop growing when (i) their effect on typical 
particles in the plasma becomes non-perturbative, as in the abelian case, or (ii) when their 
interactions with each other become non-perturbative. The first case occurs when gauge 
fields A and magnetic fields B are of order [26] 

A~P]^ an d B~ hcldP ^\ (1.1) 
9 9 

where p part is the typical momentum of the particles, and ka e \d is the wavenumber associated 
with the Weibel instability. The second case corresponds to 

and B^MsM. (1.2) 
9 9 

For weakly-interacting plasmas at times late enough that the particles have diluted to per- 
turbative densities (number densities small compared to P^x,/ a s)-> one finds that 

Afield < Ppart, (1-3) 

and so the second scale (1.2) is parametrically smaller than the first (1.1). 

Based on arguments about the form of the magnetic potential energy in anisotropic 
plasmas, Arnold and Lenaghan [26] conjectured that the fields associated with non-abelian 
plasma instabilities are dynamically driven to line up in color space at the scale (1.2) when 
their self-interactions become important, and that they then grow as approximately abelian 
configurations to the larger scale (1.1). They tested this conjecture numerically in a sim- 
plified 1+1 dimensional toy model of QCD fields in an anisotropic plasma. Subsequently, 
Rebhan, Romatschke, and Strickland [27] simulated the full hard-loop effective theory of the 
problem in 1+1 dimensions. They indeed found unabated exponential growth beyond the 
non-abelian scale (1.2), although the abelianization of the fields was not as global as that of 
the earlier toy model simulations. 

The purpose of this paper is to investigate, through simulations, whether the full 3+1 
dimensional theory behaves similarly. We find significant differences. We will soon discuss 
the details of precisely what we simulate, but here we give a preview of our results and 
conclusions. The solid line in Fig. 2 shows the growth of magnetic energy with time for a 
representative simulation. For comparison, the dashed line shows a similar 3+1 dimensional 
abelian simulation, and the dotted line a 1+1 dimensional non-abelian simulation. In the 
lower-left shaded region of the plot, the fields are perturbative, and all the curves grow at an 
exponential rate (as predicted by linearized analysis of the instability). For the non-abelian 
curves, a change takes place in the unshaded region, which turns out to roughly correspond to 
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FIG. 2: A representative simulation of instability growth for 3+1 dimensional SU(2) gauge theory 
(solid line), 3+1 dimensional abelian gauge theory (dashed line), and 1+1 dimensional SU(2) 
gauge theory (dotted line). The latter is qualitatively similar to the results of Ref. [27]. The 
parameters used for these simulations are explained in Sec. II. For 3+1 dimensions, they are 
I mH = 24, lattice spacing a = 0.25m" 1 , time step 0.1a, volume L 3 = (64a) 3 = (16m^) 3 , and 



initial amplitude A = 0.02 rno^ 1 * . For 1+1 dimensions, they are the same except that the length 

-1/2 



-1/2 

is L = 8192 a = 2048 to^ 1 and A = 0.014 rri^ 1 *- To simplify comparison of the curves, we have 
lined them up at early times by shifting the origin of time for the solid and dashed lines, which 
depend on details of initialization. 



the non-abelian scale (1.2) where self-interactions of the fields become important. The 1+1 
dimensional non-abelian simulation soon resumes exponential growth. The 3+1 dimensional 
non-abelian simulation, however, has a different large-time behavior. Though the magnetic 
energy continues to grow with time, it eventually becomes linear with time, rather than 
exponential, in the upper-right shaded region of Fig. 2. This can be seen more clearly in the 
linear- axis plot of Fig. 3. 

As we shall explain below, the hard-loop effective theory that we simulate treats the 
plasma particles as having arbitrarily large momentum, p part — > oo, and so the ultimate 
scale (1.1) limiting growth of magnetic fields is pushed off to infinity and irrelevant to 
the interpretation of our results. These simulations are designed for the sole purpose of 
cleanly studying what happens as one passes through the non-abelian scale (1.2). For similar 
reasons, the dependence of our results on g is determined by trivial scaling arguments: In the 
effective theory we use, all dependence on g can be absorbed by simple rescaling. The only 
assumptions are that g is small enough, and the separation (1.3) of scales significant enough, 
for the hard-loop effective theory to be valid. Given these assumptions, the numerical results 
we quote are valid for any g. 

The primary motivation of studying instabilities in non-Abelian plasmas is to help under- 
stand what happens in an expanding quark-gluon plasma produced in a heavy-ion collision. 
In this work, however, we simulate non-expanding systems, both for simplicity and to isolate 
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FIG. 3: Same as Fig. 2, except that the vertical axis is linear rather than logarithmic, and the time 
axis has been extended. 

the particular issue we are studying. Ignoring the expansion is a perfectly good approxima- 
tion in cases where the instability growth rate is large compared to the expansion rate. For 
example, in the context of the original "bottom-up" thermalization scenario of Baier et al. 
[9], this approximation would be valid [10] for the initial, exponential growth of instabilities 
at times late compared to the saturation time scale l/Q s which characterizes the initial mo- 
ments of the collision. 3 (The saturation momentum scale Q s characterizes the momenta of 
the original, non-perturbative, small-x gluons that eventually develop into the quark-gluon 
plasma in the saturation picture.) 

The rest of this paper gives details of our simulations. In Sec. II, we explain our formula- 
tion and discretization of the hard-loop effective theory and our choice of initial conditions. 
Sec. Ill gives further simulation results that aid in understanding the nature of the insta- 
bility in the non-perturbative regime. In Sec. IV, we discuss sources of systematic errors in 
our simulations and argue that the qualitative behavior shown in Fig. 3 is not a simulation 
artifact. Finally, we offer some last thoughts in Sec. V. Some technical results used in the 
paper are left for appendices. 



3 In this work, we will not attempt to deduce the ultimate effects of self-consistcntly including the physics 
of instabilities in the bottom-up scenario. For a recent attempt at making progress in this direction, see 
Ref. [28]. 
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II. WHAT WE SIMULATE 



A. Equations of motion 

As in the 1+1 dimensional studies of Rebhan et al. [27], the starting point for our 3+1 
dimensional simulations will be hard-loop effective theory [29, 30]. This is equivalent to 
studying kinetic theory of particles in the plasma, coupled to soft gauge fields, in an approx- 
imation where the effect of the soft fields on particle trajectories is taken to be perturbative 
[31-34]. This description is valid when there is a separation of scales fc so ft <C p pa rt as in (1.3), 
and when the soft fields have not reached the ultimate limiting amplitude (1.1). Excitations 
with momenta of order p pa rt (e.g. the initial post-collision gluons that provide the starting 
point for the formation of the quark-gluon plasma) are grouped together and described by 
a classical phase space density f(p,x,t). The softer fields associated with the instability, 
with momenta of order /c soft , are described by classical gauge fields. Studying the instability 
by treating these fields as classical is a good approximation because the instability quickly 
drives the fields to become classically large. 

Our formulation of the continuum effective theory, presented below, is equivalent to that 
used by Rebhan et al. The theory must be discretized for simulations, and we use a different 
method for discretizing velocity space than Rebhan et al, which we describe below. 

We start with the usual kinetic theory description of a collisionless plasma in terms of 
f(x, p, t) and soft gauge fields A^x, t). Then write / = fo(p) + Sf(p, x, t) and linearize the 
theory in 5f, which corresponds to a perturbative treatment of the effect of soft fields on 
hard particles. The result is well known to have the form [31] 



where 5f is in the adjoint color representation and a is an adjoint color index. Here D 
is the adjoint-representation gauge-covariant derivative. The group factor t R is defined by 
tr(TgT|.) = t R 5 ab , where T R is the color generator for the color representation R of the 
particles, and there is an implicit sum over particle species and spins in (2.1b). We use 
( — h++) metric convention. 

For ultra-relativistic plasmas, it is possible to integrate out the dependence on \p\ by 
defining 



[(A + v-D x ) 5f] a + g(E + vxB) a - V p f = 





(2.1b) 




(2.2) 



where v = p is a unit vector. Eqs. (2.1) then imply 



(A + v ■ D x ) W + gH R (E + v x B) ■ [ 



oc 




V p / = 0, 



(2.3) 



(2tt) 3 




(2.4) 



where J indicates integration over the unit sphere, normalized so that 




(2.5) 
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It turns out that these equations only depend on f through the angular function 

^'^fwT 1 ' (2 ' 6) 

which was introduced in Ref. [10]. Specifically, we show in Appendix A that (2.3) can be 
rewritten as 

(D t + v D x ) W + [E ■ (V« - 2v) - B ■ (v x V V )]M = 0, (2.7) 

where V„ is the gradient operator for the two-dimensional curved space S 2 of u's. More 
concretely V v can be related to the ordinary three-dimensional gradient by writing 

Vj, = ((J ij '-uV)|p|Vj (2.8) 

and 

V„A<(tf) = |p|V p A<Q^ • ( 2 - 9 ) 

Everything relevant about the initial distribution is specified by the angular function 
M(v). This can be split into a single dimensionful scale 

rr& = f M (2.10) 

J V 

and a dimensionless angular function 

'"oo 

The mass moo turns out to be the effective mass in the dispersion relation u> 2 ~ + m 2 ^, 
for large-momentum transverse plasmons (p ^> moo) [35]. In the isotropic case, Q = 1 and 
moo = m D I V%, where m D is the Debye mass. 

In order to discretize the problem for simulation, we generalize the procedure used by 
Bodeker, Moore, and Rummukainen [36] for the isotropic case. Specifically, we will discretize 
the 2-dimensional velocity space by expanding functions of v in spherical harmonics Yi m (v) 
and truncating at some maximum value Z max of I. We will have to check later, of course, 
that our simulation results are insensitive to the exact value of / max used, provided it is large 
enough. So, we write 

W a (v)= ]T Y. W ?™YUvl (2.12a) 

where we have chosen to normalize spherical harmonics without the usual factor of VAn, so 
that Yoo = 1 and 

Y lm (v)* Y Vm ,(v) = 8 W b mm i . (2.13) 



The bar over the Y serves as a reminder of this normalization. For working in Zm-space, 
it is convenient to rewrite the equation (2.7) for W in the equivalent form (discussed in 
Appendix A) 

(D t + v-D x )W+(E-(±[v,L 2 ]-v)-iB-L)M = (2.14) 
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where 



L= - 



iv x V v = - 



ip x V p 



(2.15) 



is the operator, analogous to angular momentum, for which the eigenvalues of L 2 and L z , 
acting on Yi m (v), are /(/ + 1) and m. One may then find the coupled equations for the Wj m 's 
using manipulations familiar from the quantum mechanics of spin. 

In this paper, we will only simulate initial particle distributions fo(p) which are axially- 
symmetric about the z axis, for which there are many simplifications. In this case, 



The equations for the Wj m 's in this case are given explicitly in Appendix A. They are the 
same as those used by Bodeker et al. [36] except for the E and B terms in the W equation 
(2.3), which are sensitive to anisotropy in f . With this exception, our discretization of the 
problem, and the evolution algorithm we use, are identical to Ref. [36]. For computational 
simplicity, we simulate SU(2) gauge theory instead of SU(3) gauge theory. We do not know 
of any reason that SU(3) would be qualitatively different. For future reference, we note that 
the vector current in Maxwell's equation (2.4) is determined by the 1=1 components of W 
as 



We should mention that, in the traditional plasma physics literature, the 1+1 dimensional 
simulations of Rebhan et al. would be referred to as 1D+3V simulations, indicating that they 
treat gauge fields and particle distributions as depending on only four of the six dimensions of 
phase space: one dimension of x space and three dimensions of momentum (velocity) space. 
All three spatial components (A x , A y , A z ) of 3-dimensional gauge fields A are simulated 
in 1D+3V, but they depend on only one spatial dimension, e.g. Ai = Ai(t,z). The 3+1 
dimensional simulations in this paper are correspondingly referred to as 3D+3V. However, 
in the ultra-relativistic limit, calculations are simplified by the fact that velocity space is 
effectively 2 dimensional, since \v\ = 1. 

Other 1D+3V simulations of non-abelian plasma instabilities have been performed by 
Dumitru and Nara [37]. Instead of working with a phase space distribution /, they simulate 
a finite number of discrete, classical particles with classical color charges interacting with 
the soft fields [38-41]. If one linearizes in the perturbations to the hard particles, this 
formulation leads to an effective theory equivalent to that above. However, they do no such 
linearization, since their interest lay in studying what happens if the effects on hard particles 
eventually become substantial. 4 



4 It is quite interesting to study the non-linearized theory, but a few caveats of interpretation should be 
kept in mind. At some time, energy loss of the hard particles through hard Brcmsstrahlung (catalyzed by 




(2.16) 



i<i 




(2.17b) 



(2.17a) 



(2.17c) 
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B. Choice of /o 



If we measure all quantities in units determined by the single dimensionful scale m^, then 
the only sensitivity of our problem to the choice of the initial hard particle distribution fo(p) 
is through the angular function Q(v). In the axi-symmetric case, that means we have to 
choose 0,(9), where 9 is the angle of v with the z axis. In this paper, we will only investigate 
results for a single choice of Q(9). Our criteria for choosing this distribution were that (i) 
perturbatively, the dominant instability should be a Weibel instability; 5 (ii) Qi m should be 
dominated by relatively low Vs, to aid in numerical convergence of our simulations to the 
^max— >oo limit; (iii) the Weibel instability growth rate should be reasonably large (relative 
to wavelength), to help minimize simulation time; (iv) Q(9) should be everywhere non- 
negative, since the distribution / is; 6 and (v) preferably Q(9) should be monotonic for 
< 9 < 7r/2. The last condition is simply superstition: We did not want to have to worry 
whether a multiple-hump distribution might perhaps have qualitatively different behavior 
than a single-hump one, and so we chose to study the simplest, most natural case. After 
some experimentation with distributions containing only / < 6 harmonics, we settled on 

n(6) oc (cos 2 /3-cos 2 ^) 3 + sin 6 /5 with (3 = 0.480. (2.18a) 

This corresponds to 

Q 00 = 1 , fi 20 = -0.790 , tt 40 = 0.367 , fi 60 = -0.093 , (2.18b) 
with all other Q tm = 0. Fig. 4 shows a plot of the angular dependence. 

2.5 

2 

1.5 

1 

0.5 

n_ n_ 3 7T 71 

4 2 

FIG. 4: The angular dependence 0.(6) of our hard-particle distribution (2.18). 

The dominant instability of this distribution is a Weibel instability with wave-number k 
along the z axis. A perturbative analysis of this instability (see Sec. IV A) shows that the 
mode with the maximum growth rate 7, in the infinite- volume, continuum limit, has 

k = 0.827moo, ^ = 0.273771^. (2.19) 

interaction with the soft fields) may become an important process. Processes which change the number 
of particles with momenta of order p par t cannot be described with the collisionlcss Boltzmann equation. 
Further, the lattice implementation via "particle and cell" codes, as used in [37, 41], suffers from spurious 
interactions between particle degrees of freedom and the most ultraviolet lattice modes, which may be 
problematic in 3D+3V simulations. For a discussion, see Ref. [41]. 

5 For a discussion, in the QCD literature, of other possibilities such as electric charge separation (Buneman) 
instabilities, see Refs. [10, 23], as well as references to the traditional plasma literature in Ref. [10]. 

6 Though this is a sensible physical requirement, it is not evident that it makes any qualitative difference 
in simulations of the harddoop effective theory. 
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C. Initial conditions 



We want to choose initial conditions that (i) have a random mix of perturbatively unstable 
modes, and (ii) are insensitive to the ultraviolet cutoff so that, for instance, the energy- 
density has a good continuum limit. A simple choice that works in any dimension is to 
choose Gaussian noise for the gauge fields A with an exponential fall-off in k: 

/d d k 
_A(fc)e*-, (2.20) 

where d—1 or 3 is the number of spatial dimensions and the A?(k) are Gaussian random 
variables with variance 7 

{A1{k)* A){k')) = (^ e ^ Vfco2 ) S ab 5 ij {27r) d 5 {d \k-k'). (2.21) 

Here A and ko are constants. It is conventional in simulations of SU(2) gauge theory to 
rescale the definition of fields to absorb factors of 2/g, but in this paper we will show 
all factors explicitly. 8 Our simulations are carried out in A =0 gauge, but the evolution 
equations and all the observables we report are gauge invariant. We choose 

E = —A = (2.22) 

and 

W = (2.23) 

as our remaining initial conditions. Both are motivated solely by simplicity. In particular, 
(2.22) automatically implements Gauss' Law. 

For small A (so that perturbation theory applies) in three dimensions, these initial con- 
ditions correspond to to an initial magnetic energy density of 

Ivk 5 A 2 

- 7 (2 - 24) 



in the infinite- volume, continuum limit, where v = 6 is the number of gauge boson degrees 
of freedom [2 spin times 3 color for SU(2)]. In 1D+3V dimensions, the corresponding result 
is 

vk\ A 2 
16(2vr)V2 g 2 ■ 

In the 3D+3V dimensional simulations reported in this paper, we choose the momentum 
cut-off scale ko in (2.21) as 

ko = 2moo. (2.26) 



^-TTtStt^- (2-25) 



7 On the lattice, we replace k 2 by ^ sin 2 (fc i a/2) in (2.21) and define the initialization of link matrices 
by Ui(x) = Gip[igaA^(x)a a /2}. 

8 That is, our normalization convention will be the traditional, perturbative convention that = — 
igA®a a /2 in the fundamental representation, where the a a are the Pauli matrices. The redefinition 

-> 2A^/g would instead make = - iA^a a . 
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Unless otherwise stated, we take the initial amplitude A in 3D+3V simulations to be 

A = 0.02m^ 1/2 . (2.27) 

In 1D+3V simulations, we take A = 0.014 m^ 2 , which corresponds to roughly the same 
value of the dimensionless 9 energy g 2 (jB 2 ) /m^. 

To perform abelian simulations, we use our SU(2) simulation code with initial conditions 
that lie in a single direction in adjoint color space, i.e. 

W(fe )4(*0)=(fe-«) 2 « % , 2 ^^( fc - fc ',. ,,28) 



III. ADDITIONAL RESULTS 

In addition to the magnetic energy plotted in Figs. 2 and 3, it is interesting to see a 
breakdown of other components of the energy of the soft fields. For an isotropic distribution 
/o, the combination 

(3.1) 



\E 2 + \B 2 + \m-J [ W 2 

J V 



would be conserved [42]. We will therefore loosely refer to ^m^ 2 f v W 2 as the U W field 
energy density." Figs. 5 and 6 compare various components of the volume-averaged energy 
density as a function of time, including the magnetic, electric, and W field energy densities. 
As one indication of the anisotropy of the soft fields, we also show \B 2 and \E 2 . Note that, 
unlike Fig. 2, we now show time all the way back to the initial conditions at t — 0. There is 
a very early transient at moot < 1 that is difficult to see in the plot, when the initial energy 
in the magnetic fields is quickly shared with the other degrees of freedom, E and W. The 
unstable modes start to grow, but it is only when they grow large enough to dominate the 
energy density that this manifests as the start of exponential growth in the magnetic energy, 
around m^t ~ 10 to 15. 

Fig. 7 shows the ratio of the B\ contribution to the magnetic energy density to the total 
B 2 . If the soft fields were isotopic, this ratio would be 1/3. This is indeed the value at 
the earliest times, due to our isotropic soft-field initial conditions. For 10 < m^t < 20, 
as the Weibel instability first starts to grow and dominate the magnetic field, B z drops 
dramatically compared to B. This is because the Weibel instability is dominated by modes 
with wavenumber k along the z axis, and magnetic fields are perpendicular to k. In the 
abelian case, that is the end of the story: B z continues to become more and more insignificant 
as the unstable modes grow. In the non-abelian case, B z = (V x A) z — ig[A x , A y ], and the 
non-abelian commutator contributes even when k is in the z direction. This contribution 
to B z will grow with time, as the gauge fields A x and A y grow due to the instability, and 
B z should become the same order of magnitude as B x and B y once those fields become 
non-perturbatively large (1.2). This growth corresponds to roughly 23 < m^t < 30 in the 
figure. At moot ~ 30, the ratio is no longer small, and so one would expect this to be where 



In d spatial dimensions, the coupling g has mass dimension (3 — d)/2 and gB has mass dimension 2. Since 
g can be scaled out of the classical equations of motion by A — > A/g, the natural dimensionless measure 
of energy for classical simulations is g 2 {\B 2 )/m^ xl rather than, for instance, {^B 2 )/m < ^ 1 . 
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FIG. 5: Various components of the energy density for the same three dimensional non-abelian 
simulation as in Fig. 2. 
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FIG. 6: Same as Fig. 5, but with a linear vertical axis. 



some deviation from the behavior of abelian instabilities should occur. This was indeed the 
case in the non-abelian magnetic energy curve of Figs. 2 and 5, where one sees a first small 
bump in the plot of magnetic energy vs. time at moot ~ 30. 

Based on arguments concerning 1+1 dimensional configurations of fields, it was conjec- 
tured in Ref. [26] that dynamics at the non-abelian scale would cause the field configura- 
tions to approximately abelianize, and that the Weibel instability would then again take 
over, causing the field to grow all the way to the ultimate scale (1.1). If abelian-like Weibel 
modes come to dominate, then B z should again become small relative to B. This was seen 
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FIG. 7: The ratio of the B 2 contribution to the magnetic energy density to the total B 2 for the 
3+1 and 1+1 dimensional simulations of Fig. 2. To make the curves line up, the time origin of the 
1+1 dimensional simulation has been shifted just as in Fig. 2 (moot — > ra^i + 9.2). 

in numerical simulations by Rebhan et al. and is reproduced in our own 1+1 dimensional 
data in Fig. 7 for moot > 30. For a brief time, the 3 + 1 dimensional simulations behave 
similarly, but the decrease of B z /B eventually stops and reverses at moot ~ 36. The ratio 
then starts growing again and finally levels out near 0.29 (slightly lower than one third) at 
moot ~ 53, which is near the beginning of the linear growth of magnetic energy in Figs. 3 
and 6. 
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FIG. 8: The ratio of electric to magnetic energy. 
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In Fig. 8, we plot the ratio of electric energy to magnetic energy. The behavior is some- 
what similar to the previous plot. Looking at the abelian curve (dashed), we see that, as the 
Weibel instability grows and dominates the energy density, the ratio drops and eventually 
asymptotes to a constant. Theoretically, this constant should be (7/A;) 2 , where 7 and k 
are the growth rate and wave number of the dominant unstable mode. Using (2.19), the 
infinite- volume, continuum value is (7/A;) 2 = 0.11, which is in good agreement with the fig- 
ure. For the non-abelian simulations, however, the ratio subsequently grows as non-abelian 
interactions first become important, and then decreases again, consistent with the picture of 
approximate abelianization. For 3+1 dimensional non-abelian simulations, however, the ra- 
tio then increases yet again, and finally approaches unity in the linear energy growth regime. 
Note that plasma oscillations with momenta large compared to the mass m^o would natu- 
rally have E 2 ~B 2 . Fig. 8 might therefore be a hint about what sort of soft field excitations 
dominate the energy in the linear regime. 

It is interesting to directly address the conjecture of Ref. [26] that non-perturbative 
dynamics would abelianize the gauge field configurations. We consider the following local 
observable, related to ones used by Refs. [26, 27]: 



c = 3 



J^((b„j,]) 2 + (b„jJ) 2 + (b' 2) jJ) 2N ' 



V2 f&\j\ 



(3.2) 



where [j x ,j y ] 2 = e abc 'j h x jy e amn -j™ j™ ', etc. The normalization of C has been chosen so that 
C would be unity if the components of j were independent random numbers with the 
same distribution. For an abelian configuration, C would be zero. Fig. 9 shows the time 
development of C in our canonical non-abelian simulations. We see that C drops suddenly 
when non-abelian interactions first become important, in agreement with the abelianization 
conjecture. However, in the 3+1 dimensional simulations, C later rises again all the way to 
unity, showing no local abelianization in the linear growth regime. 

In Figs. 7-9, we have not displayed any data for the 1 + 1 non-abelian or 3+1 abelian 
simulations at late times when the corresponding energy curves in Fig. 2 grew so large 
that they are approaching the top of the plot. At those and later times, the fields in our 
simulations grow so large that they become sensitive to the discretization of the lattice, 
and the behavior of the results is then a lattice artifact. We shall discuss this issue more 
thoroughly in Sec. IV C. 

In 1+1 dimensional simulations, the size of particle currents j was used by Rebhan et al. 
[27] to track the growth of the instability. In Figs. 10 and 11, we track a related quantity for 
our 3+1 dimensional simulations. These figures compare the magnetic energy density to the 
energy density in the 1=1 components of W. The latter is given by (3.1) as W^ m W^ m / (4771^,), 
which is directly proportional to the squared current \j\ 2 by (2.17): 

! ^w; m w; m = ^-. (3.3) 



Am 2 ^ lm lm 4m 2 



00 



Fig. 11 shows that | \j\ 2 /rn 2 ^ exceeds the magnetic energy \B 2 by roughly a factor of two 
in the exponential growth regime, but then drops substantially and is only about half the 
magnetic energy in the linear regime. 
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FIG. 9: The local measure (3.2) of the relative size of commutators, plotted as a function of time, 
for the non-abelian simulations of Fig. 2. 
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FIG. 10: A comparison of magnetic energy density (solid) and the average squared current \j\ 2 
(dashed), with the latter normalized as in (3.3), for the 3+1 dimensional non-Abelian simulation 
of Fig. 2. 
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IV. CHECKS AND SYSTEMATIC ERRORS 
A. Pert ur bat ive growth rate 

One simple check of simulations is to see whether the rate of exponential energy growth 
in the perturbative regime is consistent with the infinite-volume, continuum prediction for 
the Weibel instability growth rate. To compute the perturbative growth rate, one looks for 
exponentially growing solutions (Imu > 0) to the linear dispersion relation 

l(-u 2 + k 2 )g^ - K»K V + W(u, k)\A v = 0, (4.1) 

where K = (cu,k) and II is the hard loop self-energy [13, 22] 

dfo(p) 



IP>,fc) = g 2 t R I 
J p 



dp 1 



' v 9 + T u 

— uj + v • k — xe 



(4.2) 



In Appendix B, we formulate this in terms of the angular distribution Q(v) and then find 
a result that is well suited for calculations of growth rates in cases where f2 is expressed in 
terms of spherical harmonics. Specifically, in situations where the hard particle distribution 
fo(p) is axi-symmetric, the transverse self-energy for wave- vectors along the z axis can be 
expressed as 

U ± (u, ke z ) = v 7 ^! «, (^) (4.3a) 

with 

Kl (r,) = (1 + v 2 )$io + (1 - V 2 )[(l + l)Qi + M -(I- l)vQi(v)]- (4-3b) 
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Here, Qi(rj) is the Legendre function of the second kind defined so that it is regular at 77 = 00 
and the cut is chosen to run from — 1 to +1. 10 In cases where the dominant instability is a 
Weibel instability with k along the axis of symmetry, the dispersion relation then reduces 
to 

-u 2 + k 2 + U ± (tJ, ke z ) = 0. (4.4) 

For a given distribution Q(6), one can solve this equation numerically for each k, and then 
scan over k to look for the mode with the largest growth rate 7 = Imw. For the distribution 
(2.18) used in our simulations, the resulting 7 and k were given earlier in Eq. (2.19). 




m t 



FIG. 12: The perturbative growth in Fig. 2 compared to the predicted growth (4.5) of the dominant 
Weibel mode (2.19). As in previous figures, we have shifted the time origin of the 1+1 dimensional 
simulation to line up the curves (moo — ► moot + 9.2). 

The rate 7 gives the growth rate of the perturbative vector potentials. The corresponding 
magnetic energy should grow as the square of the field strength, so that 

\B 2 oc e 2 ^. (4.5) 

In Fig. 12, we again show our canonical simulations of Fig. 2, this time focusing on the 
perturbative regime and comparing the slope of ln(|l? 2 ) to 27. We get reasonable agreement. 
Keep in mind that there is an entire spectrum of unstable modes, not just the dominant 
mode discussed above; so the exponential growth of instabilities is not described by a single 
exponential at early times. 

B. Finite volume errors 

In Figs. 13 and 14, we show how our 3+1 dimensional non-abelian simulation results 
of Fig. 2 change if we decrease the physical volume. When comparing results for different 



For example, Q (z) = §lnf±l, Qx(z) = §lnf±i - 1, and Q 2 (z) = ^11 m £+1 _ 
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volumes, it is important to realize that different simulations will have different random initial 
conditions — there is no good way to start two simulations with the same initial conditions 
when they have different physical volumes. To get an idea of the size of this effect, we 
have simulated various volumes more than once, changing the seed of our random number 
generator. For very large volumes, one would expect the spread of the results to decrease 
with increasing volume because the volume average taken in computing the energy density 
averages over multiple spatial regions with different initial conditions. 



10 W 
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FIG. 13: Magnetic energy vs. time as in Fig. 2, but showing results for several different physical 
volumes. Multiple lines for a single volume correspond to different instantiations of the random 
initial conditions. 




FIG. 14: Same as Fig. 13, but with a linear vertical axis. 
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As one can see, the simulations shown in Fig. 14 still show a significant variation in 
results with respect to initial conditions at a volume of L 3 = (16/moo) 3 - However, the 
rates of growth in the perturbative, exponential regime and the non-perturbative, linear 
regime have clearly reached their large volume limits. Specifically, Fig. 13 shows that the 
perturbative growth rate (the slope of the curve for 10 < < 20) is not significantly 
affected in increasing volume from (8/moo) 3 to (16/moo) 3 . And Fig. 14 similarly shows 
little effect on the slope of the late-time linear growth behavior. This implies that the 
phenomenom of linear growth is not a finite- volume artifact. 

C. Finite lattice spacing errors 

1. 3+1 dimensions 

Fig. 15 shows the lattice spacing dependence of our simulations for a lattice volume of 
L 3 = (8/moo) 3 . We chose a smaller volume than that used in our canonical simulation 
of Fig. 2 so that we could push to smaller lattice spacing with our available computer 
resources. Recall from Fig. 13 that this smaller volume is adequate for reproducing the 
linear growth rate. Here and in all our simulations, the time step used in evolving the 
system is 5t = 0.1 a. 11 

10 2 | — 1 — I — 1 — I — 1 — I — 1 — I — 1 — I — 1 — i 
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FIG. 15: Magnetic energy vs. time as in Fig. 2, but showing results for several different choices of 
lattice spacing for a volume L 3 = (8/moo) 3 - The a=0.25m^ and a=0.125 m^ 1 curves are difficult 
to distinguish on this plot. 

Since we have already seen that different initial conditions can produce substantially 
different behaviors at this volume, especially at the transition between exponential and linear 



Reducing a therefore also reduces the time step. In addition, for one of our simulations, we checked that 
holding a fixed and reducing St to 0.05 a makes little (< 5%) difference to the result. In particular, the 
continued growth of energy in the linear regime does not appear to be a time discretization artifact. 
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m^t 

FIG. 16: Same as Fig. 15, but with a linear vertical axis. 



growth, it is essential in comparing lattice spacings to ensure that the initial conditions are as 
nearly identical as possible. In Figures 15 and 16, we have done this by drawing the smallest 
lattice spacing configuration randomly as described in subsection II C, and converting it into 
larger lattice spacing configurations via blocking. Besides the lattice spacing and volume, 
we otherwise use the same "canonical" values for variables as in most previous simulations. 

From Fig. 15, one sees negligible spacing dependence in the early, perturbative growth 
of the instability. From Fig. 16, we conclude that spacing does not have a significant effect 
effect on the slope of the linear growth in magnetic energy provided a is at least as small 
as our canonical value of 0.25 m^ 1 . The largest lattice spacing shows approximately the 
same exponential growth behavior, but different linear growth behavior. This is apparently 
because there is more energy in ultraviolet degrees of freedom during the linear growth 
period, and these degrees of freedom are compromised by the large lattice spacing. 

To test dependence on both lattice spacing and initial amplitude A, we show in Fig. 17 
what happens if we start with large, non-perturbative initial conditions 12 of A = 2.0. Here 
the system goes directly into linear growth, and the a = 0.25 and a = 0.125 m^ 1 results 
are reasonably close. 

From the consistency of the slope of the linear regime for a < 0.25 moo, we conclude that 
the existence of the linear growth regime is not an artifact of finite lattice-spacing. 

2. 1+1 dimensions 

To sharpen our conclusions that the linear growth regime in 3+1 dimensions is not an 
artifact of our simulations, it is useful to look at a different case, where the end of exponential 



12 For historical reasons concerning the development of our methods, this data was produced with a slightly 
different procedure for matching initial conditions between different lattice spacings: initial conditions for 
coarser lattices were obtained from those for finer lattices by setting Fourier amplitudes Aj~ exactly the 
same for the physical momenta k that exist on the coarser lattice (choosing k so that — 7r/a < k < n/a). 
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growth is a lattice artifact. Such an example is provided by our 1+1 dimensional simulations, 
shown in Fig. 18 for a variety of lattice spacings. 13 We display results to much larger 
energies than in Fig. 2. Superficially, each individual curve looks vaguely similar to our 3+1 
dimensional results: each shows an eventual end to exponential growth, although at much 
higher energy density than in 3+1 dimensions. 

In our simulations, we implement 1+1 dimensions by using our 3+1 dimensional code on a 
periodic lattice that is a single lattice spacing wide in the x and y directions. No matter how 
fine the lattice is compared to the wavelength of the unstable modes, the continuum limit 
will break down when the fields become so large that the magnetic energy per plaquette is 
of order its maximum possible value on the lattice, corresponding to ^B 2 ~ l/(g 2 a 4 ). Even 
at lower fields, the dynamics is modified by "irrelevant" operators induced by the lattice, 
such as a 2n B n+1 , which become more and more important as B grows larger. As one takes 
the lattice spacing a smaller and smaller, the fields should be able to grow larger and larger 
before these problems arise. The termination of exponential growth in Fig. 18 clearly shows 
this behavior, demonstrating it is a lattice artifact. 

In discussing our 3+1 dimensional simulations, we have emphasized how the slope of the 
late-time linear growth is not significantly sensitive to decreasing the lattice spacing. In 
contrast, a similar look at the slope of the phony late-time behavior of the 1+1 simulations, 
in Fig. 19, shows no such spacing independence. 

Fig. 20 shows the measure C of local commutators, defined by Eq. (3.2), for the 1+1 
dimensional simulations of Fig. 18. The reader may now see late-time behavior that we did 
not display in Fig. 9, and also see that this behavior is a lattice spacing artifact. This is the 
reason we truncated our 1+1 dimensional curves in Figs. 7-9. Specifically, we truncated 
those curves when the magnetic energy densities reached 20 m^/g 2 , which is a little before 
where the a = 0.25 m^ 1 energy curve deviates from the continuum limit in Fig. 18 and 
where the corresponding curve for C suddenly begins to rise in Fig. 9. A similar truncation 
was made for the Abelian results, where a similar large-field issue arises since we implement 
compact rather than non-compact Abelian gauge theory. 



Again for historical reasons, initial conditions were matched using the method of footnote 12. 
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FIG. 18: Magnetic energy vs. time for 1+1 dimensional simulations with several different choices 
of lattice spacing for a system length L = 128/moo. The remaining parameters are as in Fig. 2. 
except that here we have not shifted the origin of the time axis. 
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FIG. 19: Same as Fig. 18, but the vertical axis is linear and different vertical offsets have been 
added to each curve to fit them onto the plot. From top to bottom on the left side of the plot, the 
curves correspond to arrioo = 0.125, 0.25, 0.5, 0.0625, and 0.03125. 



D. Finite Z max errors 



Figs. 21 and 22 show the dependence of our results on Z max . The figures show simulations 
with fixed lattice spacing a = 0.25 and fixed physical volume L 3 = (8/moo) 3 . By 
comparison of the Z max =24 and Z max =48 curves in Fig. 22, we conclude that our canonical 
choice of Z max = 24 seems to be a good approximation to the large Z max limit and that the 
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FIG. 20: The local measure (3.2) of the relative size of commutators for the various 1+1 dimensional 
simulations of Fig. 18. From top to bottom, the curves are in descending order in the size of a. 




FIG. 21: Magnetic energy vs. time for several different choices of / max . These simulations have 
a = 0.25m" 1 , Z max = 24, L 3 = (32a) 3 = {8m^f, and A = 0.02 m^ /2 . The Z max =24 and / max =48 
curves are difficult to distinguish on this plot. 



linear growth regime is not a finite Z max artifact. 

We have chosen to run simulations only for even values of / max , based on the experience 
of Ref. [36] for isotropic particle distributions, where it was found that convergence to the 
large Z max limit was much faster for Z max even than / max odd. 
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FIG. 22: Same as Fig. 21, but with a linear vertical axis. The I 
again very close together. 



max — 24 and l ma x=48 curves are 



V. DISCUSSION 



We have found that 3D+3V dimensional non-abelian plasma instabilities behave qualita- 
tively differently than 1D+3V dimensional instabilities once they grow non-perturbatively 
large. Initially, there is a period of non-perturbative growth that looks quite similar to the 
1D+3V case and to early conjectures about abelianization, but eventually the 3D+3V in- 
stabilities settle into a period characterized by linear rather than exponential growth of the 
magnetic energy. Based on our analysis of sources of systematic error, we believe that this 
conclusion is not a simulation artifact. 

There are many more things one would like to know about the linear growth regime, such 
as its efficiency at scattering and isotropizing typical particles in the plasma, its power spec- 
trum, and possible models for the underlying physical processes. It would also be useful to 
check whether non-perturbative linear growth occurs for a wide variety of different distribu- 
tions /o(0), beyond the single case studied here. We leave this and further characterization 
to future work. 
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APPENDIX A: HARD LOOP EFFECTIVE THEORY IN im-SPACE 

First, we will fill in a few steps in the text. The transition from (2.3) to (2.7) follows by 
using (2.8) to write 

V p /o = -V„/ + ^. (Al) 
V op 



24 



Then integrate the last term by parts in dp to obtain 

and so (2.7). Using (2.8) and (2.15), one may derive that 

[v,L 2 } = 2(V v -v), 



f 

Jo 



(A2) 



(A3) 



which yields (2.14). 

Now use the /m-expansions (2.12) of W and Q in (2.14) and project out the /m-th 
component of the result. Working in A =0 gauge, which is the gauge of our simulations, 
this yields the evolution equation for Wi m , 

d t W lm = J2{-( lm \v\l'm') ■ D x W llm/ 



(E ■ (lm\\[v, L 2 } - v\l'm') - iB ■ (lm\L\l'm')) m^JW} 



= ^2{-(lm\v\l'rri) ■ D x W Vm , + 



1 + 



1(1+1) -l'(l'+l) 



E ■ (lm\v\l'rri) tt{ 



+ im 2 ^ B ■ (lm\L\l'm')Qii m i^ . 



(A4) 



The expectation values (lm\v\l'm') and (lm\L\l'mf) are simple results from the quantum 
mechanics of spin. The first, (lm\v\l'm') , is also relevant to the simulations for isotropic /o, 
and explicit formulas may be found in Appendix A of Ref . [36] . The L expectations are 

i(lm\L x \l'm') = ^5 tt > 5 m ,m'+i \/(l + m)(l- m 1 ) + 5 m > jTn+1 \/(l + rri) (I - m) , (A5) 
i (lm\L y \l'm') = 6 m ,m'+i a/(/ + m) (/ - m') - S m > jm+1 a/(/ + m!) (I - m) , (A6) 



i (lm\L z \l'mf) = im5u'5 ri 



(A7) 



In practice, we found it convenient to implement (A4) with a basis of real functions of v 
instead of the usual complex F; m 's. So we switched to the basis of 



V2ReY lm , m>0; 

Yim = { Y tm , m = 0; 

y/2 1mYi\ m \, m < 0, 

with overall normalization again set by (2.13). 



(A8) 



APPENDIX B: SELF ENERGY II 



In this appendix, we derive the result (4.3) for the hard-loop transverse gluon self-energy 
n_i_(u;, k) when the particle distribution is axi-symmetric and k is along the axis of symmetry. 
By integration by parts, the spatial part of Eq. (4.2) can be recast into the form [10]: 



n ij (u,k) =e 2 



Hp) 



p 



9' 



k i v j + k j v i 



+ 



(-co 2 + Ar>V 



— oj + v ■ k — ie (-co + v ■ k — ie)' 



(Bl) 
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Take the axis of symmetry to be the z axis. For k along that axis, axial symmetry implies 

n_L = Hxx = n ra = — kikj)iiij. (B2) 

Then, using (Bl) and introducing r\ = cu/k, 



ilj_(a;,A;e 2 ) = e z 



i: 



V 



d(cos 6) 



1 + 



1 + 



;i-^)[i-Kfc) 2 ] 

2(-r] + v- k) 2 

1 - T] 2 )[l - COS 2 e\ 



2(-r] + cos6f 



M{0) 



Expanding Ai in spherical harmonics gives 



Pi(x) 



Decomposing by partial fractions and using 



J\x 



Pi(x) 

-1] + X 



and its 77- derivative 



(-77 + x) 2 
produces the final result (4.3). 



d, Pl{x \„ = -205(7,) = 2 -p^- [Qi + M -vQiiv)} 



l-T] 2 



(B3) 



(B4) 



(B5) 
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